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Abstract 

Inelastic scattering and carrier capture by defects in semiconductors are the primary causes of hot- 
electron-mediated degradation of power devices, which holds up their commercial development. At the 
same time, carrier capture is a major issue in the performance of solar cells and light-emitting diodes. A 
theory of nonradiative (multiphonon) inelastic scattering by defects, however, is non-existent, while the the¬ 
ory for carrier capture by defects has had a long and arduous history. Here we report the construction of a 
comprehensive theory of inelastic scattering by defects, with carrier capture being a special case. We distin¬ 
guish between capture under thermal equilibrium conditions and capture under non-equilibrium conditions, 
e.g., in the presence of electrical current or hot carriers where carriers undergo scattering by defects and 
are described by a mean free path. In the thermal-equilibrium case, capture is mediated by a non-adiabatic 
perturbation Hamiltonian, originally identified by Huang and Rhys and by Kubo, which is equal to linear 
electron-phonon coupling to first order. In the non-equilibrium case, we demonstrate that the primary cap¬ 
ture mechanism is within the Bom-Oppenheimer approximation (adiabatic transitions), with coupling to the 
defect potential inducing Franck-Condon electronic transitions, followed by multiphonon dissipation of the 
transition energy, while the non-adiabatic terms are of secondary importance (they scale with the inverse 
of the mass of typical atoms in the defect complex). We report first-principles density-functional-theory 
calculations of the capture cross section for a prototype defect using the Projector-Augmented-Wave which 
allows us to employ all-electron wavefunctions. We adopt a Monte Carlo scheme to sample multiphonon 
configurations and obtain converged results. The theory and the results represent a foundation upon which 
to build engineering-level models for hot-electron degradation of power devices and the performance of 
solar cells and light-emitting diodes. 
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I. INTRODUCTION 


Elastic scattering of electrons by phonons, impurities, and other defects limits the conductivity 
in metals and the carrier mobility in semiconductors. The fundamental theory is well established, 
parameter-free mobility calculations have become possible [[11,121], and engineering-level modeling 
methods are widely available. Inelastic scattering of hot electrons by defects has long been known 
to cause device degradation. For example, hot electrons in Si-Si02 structures can transfer energy 
and release hydrogen from passivated interfacial Si dangling bonds l^lSO. More recently, it was 
found that hot electrons cause degradation of power devices based on wide-band-gap semicon¬ 
ductors [1^. It has been shown that the degradation is caused by hot-electron-mediated release of 
hydrogen from hydrogenated defects such as Ga vacancies or impurities [TD. In other cases, carrier 
(^ture transforms benign defects to metastable configurations that cause recoverable degradation 
iSn . Similarly, non-radiative carrier capture by defects, which is a special case of inelastic scatter¬ 


ing, limits the performance of photovoltaic cells, light-emitting diodes and other devices la. 


m 


A theory of inelastic scattering by defects by multiphonon processes (MPPs) does not exist 
while the theory of non-radiative carrier capture or emission by defects by MPPs has a long and 
controversial history. In 1950, Huang and Rhys 111 ill reported a theory of how the energy of 
lattice relaxation that accompanies the photoionization of a defect is dissipated by MPPs. The 
process was described within the Bom-Oppenheimer or adiabatic approximation (BOA) and the 
Frank-Condon approximation (FCA). The former says that the electronic and nuclear (vibrational) 
wave functions obey decoupled equations. The latter states that an electronic excitation occurs 
instantaneously and relaxation processes follow at a relatively slow pace, allowing one to write the 
excitation rate (Fermi’s golden rule) as a product P = AF, where A describes the instantaneous 
electronic excitation in the initial lattice configuration and F, the so-called line-shape function, 
describes the MPPs that occur during lattice relaxation. In the Huang-Rhys theory, the operator 
that causes the excitation is strictly the photon field and MPPs dissipate only the energy of the 
ensuing lattice relaxation. 


In the same paper, Huang and Rhys llllh also proposed a theory for non-radiative multiphonon 


transitions between defect levels. Such transitions are caused by the terms that are dropped when 
the Bom-Oppenheimer approximation (BOA) is made, namely derivatives of the electronic wave- 
functions with respect to nuclear positions (non-adiabatic terms). In 1952, Kubo 1121] indepen¬ 
dently invoked the same non-adiabatic terms as being responsible for the thermal ionization of a 
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defect. In subsequent years, Kubo and Toyozawa lll3ll and later Gummel and Lax [Il4l] adopted 


Kubo’s formalism to e xp 


varskii and Sinyavskii IfS 


ore carrier capture and emission using analytical approximations. Ko- 
4l7l] published several papers expanding on Kubo’s formalism. In 1977, 


in search of a practical scheme to model electron capture in experiments , Henry and Lang [Il8h 
adopted a Huang-Rhys analog: the electronic transition is caused instantaneously by the perturba¬ 
tion potential AV generated by atomic vibrations - the linear electron-phonon coupling potential 
that is normally thought to cause elastic scattering and is used for mobility calculations. The 
following year, Ridley H] showed that the Henry-Lang model exhibits the correct temperature 
dependence at high temperatures (the semi-classical limit), but pointed out that the correct way 
to calculate non-radiative capture cross sections is through the non-adiabatic perturbation terms 
identified by Huang and Rhys 111 IH and by Kubo lll2ll . In 1981, however, Huang showed that 
the non-adiabatic perturbation Hamiltonian and the linear electron-phonon coupling perturbation 
Hamiltonian are equivalent to first order ll20ll . The issue whether such a first-order calculation is 
adequate remained open as, throughout the years of all these developments, only model calcu¬ 
lations were pursued, largely analytical, employing model defect wave functions. Furthermore, 
calculations of the line-shape function were typically restricted by the assumption that a single vi¬ 
brational mode contributes to the MPPs. In the chemical literature, noradiative transitions between 


molecular orbitals have been studied 121 , 


2211 . It was recognized that inclusion of all vibrational 


modes in the MPP calculation leads to exploding computational requirements as the size of the 
molecule increases 112ill . The so-called parallel-mode approximation or simply a single vibrational 
mode are typically used ll22ll . 

The first application of modern density-functional-theory (DFT) calculations to MPPs in the 
case of luminescence, i.e., the classic Huang-Rhys problem where an electronic transition is caused 
by the photon field and MPPs dissipate the ensuing lattice relaxation, was reported by Alkauskas 
et al. 1 I 23 II . These authors studied the luminescence spectra of defects in GaN employing DFT 


pseudo wave functions for the electronic matrix elements and the single-phonon-moc 


tion to the Huang-Rhys line-shape function. In a more recent paper, Alkauskas et al. ll24ll reported 


e aproxima- 


calculations of non-radiative capture of carriers by defects using the linear electron-phonon cou¬ 
pling perturbation Hamiltonian, pseudo wave functions, and a single-phonon-mode to calculate 
the MPPs that dissipate the transition energy. They pointed out that the electronic transition is a 
slow process because capture is mediated by the phonons that are localized around the defect. 

In this paper we first revisit the theory of carrier capture by defects. We identify two distinct 
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regimes that are governed by different processes. One is carrier capture under thermal equilibrium 
conditions, i.e., capture occurs in tandem with emission and electrons in the conduction band (or 
holes in the valence band) are not being accelerated. Under these conditions, capture and emission 
are inverse processes, i.e., the role of the initial and final states is reversed. For an electron bound 
at a defect, emission amounts to a transition to a band state that is an eigenstate of the same 
Hamiltoninan (perfect crystal plus defect potential). Band states are occupied according to the 
Fermi-Dirac distribution function. Any of these carriers can be captured into the defect’s ground 
state. Under such conditions, band carriers are effectively undergoing diffusive Brownian motion. 
In this case, the Huang-Rhys-Kubo (HRK) non-adiabatic Hamiltonian perturbation is the only 
possible cause for these thermal transitions. 

Under non-equilibrium conditions, however, e.g., in the presence of an electrical current, car¬ 
riers are accelerated in a specific direction and a mean free path is defined by scattering events. It 
is then standard procedure to treat the band electrons as being in eigenstates of the perfect crystal 
Hamiltonian and consider scattering by the defects. In particular, one considers elastic scattering 
by defects as a mechanism that limits the carrier mobility. In this case, the initial and final states 
are eigenstates of the perfect crystal Hamiltonian and the defect potential acts as the perturbation 
that causes the transitions, i.e., the defect potential is “turned on” in order to use time-dependent 
perturbation theory and arrive at Fermi’s golden rule. Clearly, hot carriers can undergo inelastic 
scattering as well, dropping to a Bloch state of lower energy, with the energy dissipated by MPR 
For such calculations, one must again “turn on” the defect potential, though the HRK non-adiabatic 
perturbation must also be included. Transitions caused by the defect potential are within the BO 
approximation, whereas those caused by the HRK perturbation Hamiltonian are non-adaiabatic. 
Finally, under such non-equilibrium conditions, carrier capture can be viewed as a special case of 
inelastic scattering: if the defect potential can cause elastic scattering and inelastic scattering with 
energy dissipation via MPP, then it certainly should also be included as a cause for capture. 

In the capture case, however, there is a subtle difficulty. In order to derive a transition rate 
using Fermi’s golden rule, initial and final states must be eigenstates of the same Hamiltonian. In 
the carrier capture case, however, the final state is an eigenstate of the crystal Hamiltonian plus 
the defect potential, whereas the initial state is an eigenstate of the perfect crystal Hamiltonian. 
The difficuty can be overcome if we prepare a propagating state for the incoming electron that is 
not aware of the bound state’s existence, with capture being triggered by the sudden turning on of 
a suitable coupling (initial and final states must belong to the same Hamiltonian for the concept 


5 



of a transition to be meaningful) to the defect potential. Such adiabatic transitions have not been 
considered so far in the context of multiphonon transitions at defects in semiconductors, but they 
are commonly invoked in chemistry for elecron transitions in molecules [l25l-l27n. 

We will develop a comprehensive theory of inelastic scattering and capture for transitions 
caused by both the defect potential (adiabatic transitions) and by the non-adiabatic HRK perturba¬ 
tion Hamiltonian. We will show that, for carrier capture, adiabatic transitions are the zeroth-order 
term in an expansion in the defect-atom displacements that following capture (lattice relaxation) 
and are, therefore, dominant under non-equilibrium conditions. The electronic transition is caused 
instantaneously by the defect potential (it is effectively a Franck-Condon transition) and the en¬ 
ergy is dissipated by MPP. The next order in the series, which is linear in the atomic displacements. 


comprises two terms, only one of which has been captured by prior theories [12(1 


2411 . We estimate 


that these “linear terms” make smaller contributions to the capture rate as they scale with 1/m, 
where m is a typical nuclear mass in the defect complex. The adiabatic perturbation Hamiltonian 
that couples the incoming electron to the defect is constructed in terms of Hamiltonian matrices as 
in the Forster theory of electron and exciton transfer in molecules 112511 . which allows the derivation 
of Fermi’s golden rule for these transitions. 

In addition to presenting the basic elements of the fundamental theory, we report explicit calcu¬ 
lations for capture cross sections as functions of energy transfer for a prototype defect using DFT 
for the electronic matrix elements. We employ the Projector-Augmented Wave (PAW) scheme 


Il28ll . which allows the use of the all-electron defect potential and wave functions as opposed to 
pseudopotentials and pseudo wave functions. For the calculation of the line-shape function, we 
introduce a Monte Carlo scheme to sample the space of phonon combinations that contribute to 
the MPP energy dissipation and find that random configurations containing up to twelve different 
phonon modes and trillions of configurations are needed to obtain converged results. 

A few more observations are in order before we describe the present theory in detail. In a 
perfect crystal without defects, the HRK perturbation Hamiltonian is responsible for electron- 
phonon scattering (only linear coupling is usually included) and for the formation of polarons, 
which are electrons or holes dressed by phonons. Under strong-coupling conditions, the HRK 
Hamiltonian can be responsible for polaron self-trapping. When a defect is present, the HRK 
Hamiltonian can cause carrier capture. As Alkauskas et al. ||24ll pointed out, such capture is 
very slow. Indeed it is caused by the derivatives of the electronic wave functions with respect to 
nuclear displacements, which amounts to a “frozen electron approximation” (recall that the BO 
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approximation is effectively a “frozen nuclei approximation”). As we already noted, this kind of 
capture occurs under thermal equilibrium conditions, which corresponds to constant emission and 
capture by inverse processes, i.e., the band electrons are definitely “aware” of the defects, i.e., they 
should not be treated as “free” carriers with a mean free path, undergoing scattering by defects 
and phonons. In this regard, the linear coupling approximation [|24ll should be viewed as the zero 
mean-free-path limit, whereas the theory put forward in this paper represents the limit in which 
the mean-free-path is only bounded by ^capture, the mean distance an electron travels before being 
captured by a defect. 

The conditions under which capture cross sections are measured by junction capacitance meth¬ 


ods [IlSll are close to equilibrium, i.e., they are slow. Similarly, in light-emitting diodes, carriers 


by design have minimal acceleration through the pn-junction. However, even in such deliberate 
setups, there must still be some nonequilibrium driving forces, e.g., a current must flow through 
the system, in order to carry out the measurement or for the device to operate. The carrier mean- 
free-path is always finite, never exactly zero. Therefore, a realistic model of the measured capture 
cross sections can be obtained by scaling the difference between the two limits according to the 
factor L /Lcapture where L is the elastic scattering mean-free-path, 

L 


a = 


■ ^adiabatic T ^^nonadiabatic) 


( 1 ) 


^capture 


where Unonadiabatic is the capture cross section due to the HRK Hamiltonian and Uadiabatic is the 
adiabatic capture cross section calculated in this paper. 

For scattering of a carrier into another propagating state at a lower energy, the defect is left in 
the same charge state, which requires that scattering by the defect potential is elastic (no energy can 
be dissipated in the Franck-Condon approximation in such a case). We find that inelastic scattering 
can still occur within the BOA by the first-order correction to the Franck-Condon approximation, 
which are the linear terms discussed above. 


II. FERMI GOLDEN RULE FOR ADIABATIC AND NON-ADIABATIC TRANSITIONS 

As discussed in the previous section, in order to describe transitions, it is always necessary 
to identify the piece of the total Hamiltonian that causes the transition between eigenstates of an 
approximate Hamiltonian. Let us be more specific. In the hydrogen atom, one usually includes 
only the Coulombic attraction between the proton and the electron, leaving out the electromagnetic 
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field at large. The calculated energy levels are only eigenstates of this approximate Hamiltonian. 
The electromagnetic field, treated as a perturbation, then causes a transition from, say, a 2p state to 
the Is state. In Auger transitions, one must leave out specific electron-electron interactions that are 
then introduced to cause transitions OOl] . Our task here is to identify the approximate Hamiltonian 
whose eigenstates are the propagating state of the incoming electron that is not aware of the bound 
state of the defect potential and the final state, which can be either another propagating state that 
is not aware of the existence of a bound state at a lower energy or the bound state itself, and 
determine the perturbation Hamiltonian that causes the transition. 

In the BOA, the many-electron Hamiltonian depends parametrically on the nuclear positions 
and the total wave functions are products of many-electron wave functions and phonon wave func¬ 
tions. Within DFT, the many-electron wave functions are Slater determinants of Kohn-Sham wave 
functions. We start by defining the many-electron Hamiltonian for the perfect crystal and the 
corresponding eigenvalue problem. 


In/ nIn/ 


( 2 ) 


For the crystal containing a single defect, we have 




(3) 


One normally writes 

FT = + A/7. (4) 

The partitioning of the total Hamiltonian H according to Eq. dH) is not useful for our purposes. 
Instead, we write 

H = H^ + //f (5) 

where, 

= enl^n). ( 6 ) 

In order to obtain an explicit description of , which then defines through Eq. ([5]), we 
express AH in terms of the complete set of functions \l/„ : 

AH = J2 \^m){^m\AH | (^„ |. (7) 

m n mn 

We then define by 

= \^,)AH,f{^f\ + \^f)AHf,{^,\, ( 8 ) 
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where the subseripts i and / denote the eigenstates of that are the initial and final states of 
our problem. This definition of is analogous to the so-ealled Forster transition often used in 
energy transfer in moleeules 11251] . In effeet, eliminates the eoupling of the ineoming eleetron 
via the defect potential to the final state, whether propagating or bound. The defect potential AH, 
which can be arbitrarily strong, is still present. It is the perturbation Hamiltonian that is 
weak and can cause transitions whose rate is describable by Fermi’s gold rule, i.e., to first order in 
Note also that the state lihj) contains an incoming electron that “sees” the defect potential, 
but does not couple to the bound state. Also, for all practical purposes, for carrier capture we have 
f) = |<I)/) (i.e., the bound sate is not affected by the presence of an incoming electron that does 
not couple to the defect). 

The adiabatic transition rate is given by the usual Fermi’s golden rule by 

271 


- 

^^f - ^ 


- J2\{Xf\{^f\Hr\%m)\' SiQf -Q^ + 


(9) 


where 0jj are the total phonon energies of states \Xij) and e*/ = e/ — e* is the energy difference 
between the electronic states and |1I//). For capture, it is usually assumed that there is one 
final electronic state with a given energy difference e,/, but there are many phonon configurations 
that can make up this difference. If there are multiple electronic states at the same energy we need 
to sum Eq. @ over all such states. 

In addition to , there are terms beyond the BOA, usually referred to as the non-adiabatic 
terms |il[ 12], that cause multiphonon transitions. These terms contain derivatives of the electron 
wave functions with respect to nuclear coordinates {Rfc} and are the terms neglected when one 
invokes the BOA. They contribute to the total transition rate Wij via the matrix element. 


k 

where mk is the mass of atom k. This contribution will be discussed in detail later. 
One can define a cross section for inelastic scattering or carrier capture by 

WifQ 


( 10 ) 


( 11 ) 


where Vg is the group velocity of the incident electron, O is the volume over which the state |i) is 
normalized, so that Vg/il represents the flux of the incoming electrons. 

We will work within DFT so that the many-electron wavefunctions are Slater determinants of 
Kohn-Sham one-electron wavefunctions and the many-electron Hamiltonians are those of non- 
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interacting Kohn-Sham quasi-particles in the presence of an effective single-particle external po¬ 
tential. From now on we will view the Hamiltonians and wavefunctions in Eqs. dH) and (fT^ as 
one-electron Kohn-Sham Hamiltonians and electron wave functions without change of notation. 


A. Adiabatic series 


We now examine the electronic part of the transition matrix element in the BOA by showing 
explicitly its dependence on the atomic coordinates, 

Aff°({R,}) = |(<I>,({R,})|/ff°({R,})|'r.({R,})>p. (12) 


The BOA by itself does not separate electron and phonon matrix elements. A further approxima¬ 
tion is needed. We expand 


MfO({Rj}) = MfO({Rf}) + y^(R, - Rf) ■ VR.Aff°({R,}) + ..., (13) 

k 


in terms of the atomic displacements where R^°^ are the atomic positions in a reference 

state, which will be determined later. The transition rate is then, 

Aff°({Rf |{Af;|A',)|"3(e, - e. + £.^) 

/ 



2 

^ Vk. Aff°({Rf}). {X,|(Rj - Rf )|Afj) 3(0, - Oj + + ..(J4) 

k 

Here the cross terms are dropped because the zeroth order and first order terms cannot have the 
same final phonon wave functions - the number of phonons needed to ensure a nonzero overlap 
matrix element are different for the two cases. The first term in this expansion represents a com¬ 
plete separation of the electron and phonon wave functions as if they are independent of each other 
and corresponds to the Frank-Condon approximation. The second term is the first order correction 
to the Frank-Condon approximation arising from the BOA perturbation Hamiltonian . 



B. Non-adiabatic series 


According to Huang [|20ll . the non-adiabatic matrix element defined in Eq. (fT^ can be evaluated 
for linear phonon coupling. 


^(^/({Rf })|VR.ff.({Rf })|>I-i({Rf })) . {A',|(R, - R<“')|A'), (15) 
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where He is the electron part of the Hamiltonian. When electron-phonon coupling Hep = 
/fe({Rj}) — is introduced, the electron wave functions are changed by a perturba¬ 

tion, 

|i».({R,})) = (16) 

i'^i 


(and a similar equation for the final states). We write both the initial and final states in the form 




Substituting this into Eq (fTOl) and keeping only the linear terms. 


})l' 5 >t.>|.V)) - (.Y^|{'I'/({Rf 


(e.-e;){A';|{*;({Rf>})|i®i)|Xi) 


Gif 


(>i.,|VR.R|».)(>i-,({R'“’})l»i'({Rr})) ■ WI(R* - RDI-Y.) 


5 : 5 :, 

1 • j ■ '~'Z '~'Z 

k v^i 

. (A/KRj - r™)|a:.>. 


(0)^ 


,(0)^ 


(17) 


(18) 


Here the first equality results from the Schrodinger equations for the phonon wave functions and 
for the second equality we used 0* — 0/ = e*/. 

We note that the above linear-order term in the non-adiabatic series has the same phonon matrix 
element as the linear-order term in the BOA series of the previous section. This indicates that the 
leading non-adiabatic term is a smaller contribution to the electron capture rate compared to the 
zeroth-order BOA term. The electronic matrix element in the non-adiabatic series is different than 
the BOA series. We will show later that both these terms scale as 1/m, where m is the mass of a 
typical atom in the defect complex. 

The linear term in Eq. (fTSl) is usually referred to as the linear electron-phonon coupling term. 
A similar term has been calculated by Alkauskas et al [|24ll , with the exception that in that work 
the wave functions are | $*(/)) which are the eigenstates of the full Hamiltonian He, whereas in our 
case the wave functions are which are the eigenstates of the Hamiltonian We recover 
the term calculated by Alkauskas et al. if we combine the BOA and the non-adiabatic series. We 
make use of the result in Eq. (l2^ and get for our final result 
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h 
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Mf°({Rf })| |{A',|.Y.)|‘^i(e^ - e. + ii,) 




E [(«/|VR.-ffe|1>i> - («/l't?>(1>/|VR.ff.|4,)] ■ (AVKRi - Rf )|A'j) 


X 


X 5(0/ — ©i + eif) + ... 


(19) 


Here the first term is the zeroth-rder term that eorresponds to the Franek-Condon approximation 
and thes seeond terms is the totality of contributions from the linear terms in the two series. The 


first term in square brackets is precisely the term that Alkauskas et al. [|24ll calculated. We note 


that there exists a second term, which has the appearance of a force term. These two terms can 
either add or subtract. We will show shortly that these linear-order terms are proportional to 1/m, 
where m is a typical atomic mass in the defect complex, and are, therefore, significantly smaller 
than the zeroth-order Franek-Condon term, which is dominant. 


III. ELECTRON MATRIX ELEMENTS 

We first consider the zeroth order term in the BOA series, which yields a capture cross section 
that can be written in the familiar factorized form. 


O'if — T^ifFif, 


( 20 ) 


where A*/ contains the electronic part of the matrix element, 




(4<,({R'“>})|Rr({Rf'})|'I<i({R'“>})) 


and F is called the line shape factor due to vibrations, 

/ 

Next we will consider these two factors separately. 

Detailed derivations given in Appendices U and |n] find the final results 


( 21 ) 


( 22 ) 


= -($/|vl//)e,/, 


(23) 


and 
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VR.Aff° + (<I-,|VR.i/|<I-,) = (4/|VR.i/|4,) - |4/). (24) 

For the evaluation of the above matrix elements, we employ the PAW seheme, whieh allows us 
to use all-eleetron wave funetions instead of pseudo wave funetions. Details are given in Appendix 

nni 


IV. PHONON MATRIX ELEMENTS 


First, we eonsider the effeet of displaeements for a elassieal Hamiltonian. We derive this Hamil¬ 
tonian for the ion motion from whieh the phonon wave functions and matrix elements can be 
calculated. For this purpose we start with a supercell containing tIq number of atoms with the 
defect site at its center. This supercell is repeated N times using the Bom-von-Karman periodic 
boundary condition. For the initial state, the equilibrium positions of the atoms are Rk where the 
subscript k runs through both the atomic index within the supercell and the cartesian components. 
Each atom oscillates around its equilibrium position with displacement Uki, where the subscript 
I labels different copies of the supercell under the Born-von-Karman periodicity. Using the har¬ 
monic approximation for the potential energy, under which only terms that are second order in 
displacements make a contribution and introducing force constants we can write [|25n . 

2 




kl 


1 

—ruk 

2 


du 


kl 


dt 




2N 


k'l'Uk'V 


k'V 


(25) 


where the atomic mass ruk also carries the subscript k for convenience even though it depends 
only on the atomic index and not the coordinate component index. 

When an electron is absorbed or emitted from the lattice, the equilibrium position of the atoms 
change. The new equilibrium positions are Rk + A^. The new Hamiltonian has the same form 
after initial displacement vectors Uki are replaced by u'j^i = Uki — A^.. The final state Hamiltonian 
is then written as 

2 




kl 


djuki - Afc) 
dt 


+ ^ {uki — Afc) ^ki,k'i' {uk'v — Afc'^ 


(26) 


k'V 


where we make an assumption that force constants do not change due to the electron capture 
or absorption. Since displacements A^ do not depend on time, the kinetic energy term remains 
unchanged. Expanding the potential energy to first order in displacements reproduces the same 
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term in the original Hamiltonian plus a term that includes Uk'i'^k- 


Hf — H[ — — ^ ^ki^k'i'^kUk'i' (27) 

' kl,k'V 

Transforming to the normal-mode representation in terms of the generalized coordinates, 


= 


^ ^ \/^^k^kl'^j,klj 


(28) 


M 

where is the kith element of the eigenvector for mode j. Note that in this definition of the 
generalized coordinate qj, it has absorbed the mass factor ^/fnk. The Hamiltonian is expressed as, 


Hf — 2 Qj + 2 ^ Dkk'{^j)wjk'^/Tnk^k■, (29) 

3 3 ^ 3 kk' 

where Uj are the eigenfrequencies. A phase factor of the form exp(ikj ■ r;/), where is the wave 
vector of mode j, from Wj^k'v is absorbed into the force constant matrix $ yielding the dynamical 
matrix D, and reducing Wj^k'v to Wjk' (independent of /')• Since we assume that force constants 
remain the same after electron capture. 


Dkk'{^j)wjk> = ‘^]wjk (30) 

k’ 

The linear term causes a general coordinate displacement. 


6qj 


y/N 


^ ^ \/Hlk^kWjk- 


(31) 


We can express the normal coordinates of the lattice for the final (/) state, gj, in terms of those for 
the initial (i) state, qj, 

(If,3 = (32) 


so that the final Hamiltonian is: 


= + (33) 

3 3 

A. Zeroth-order phonon matrix elements 

We have derived the expression for the generalized coordinates resulting from the lattice dis¬ 
placements. These generalized displacements enter the phonon wave functions \X^i {qj)) and 
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\Xnf {qj + Sqj)), respectively in the quantized versions of the harmonic oscillator Hamiltonians 
HI and H'r . Now we turn to the evaluation of phonon matrix elements {X f {qj + 6qj)\X^i{qj)). 
When the displacement 6qj are small, we can show that the dominant contribution comes from 
single phonon emission or absorption for each normal mode. Suppose the initial state of mode j 
has n phonons and its final state has n + p phonons, (we dropped the index for the mode, since it 
is present in the notation of generalized coordinate). Using the integrals provided in Appendix HVl 
the matrix elements for the phonon part are, 

{A'„+i(9, + ^%)|A'„(®)) = 

+ i5®)|A:„fe)) = (35) 

The integrals for phonon modes that maintain the same occupation numbers are calculated to 
second order in qj, 

{Xn{qj + 6qj)\Xn{qj)) = 1 - (36) 

Now we consider how to evaluate Eq. (l22l) . The total number of phonon modes in the supercell 
is M = 3{na — 1) excluding the translational motion, and the total number of phonon modes in the 
entire system is MN, since supercell is repeated A^-times. We assume that there is a one-to-one 
correspondence between phonon bands before and after the capture. The wave function of the 
initial phonon state is 

MN 

|X,) = (37) 

j=i 

and that of any one of the final phonon states is 

MN 

i.Y;)=nK')' 

i=i 

where n* and nj are the occupation numbers of phonon mode j before and after the capture, 
and are also used to label the wave functions. The total phonon energies for initial and final 
configurations are 

MN 

(39) 

i=i 

and 

^ MN 

= (40) 

i=i 
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respectively, where ce* and cej is the phonon frequency of mode j in the initial and final configura¬ 
tion of the defect, respectively. With the overlap matrix for each individual mode expressed as Eq. 
([U) and using Eqs. (iJTl) . (I38l) . (1^ . and (ITO) . Eq. (I22l) now takes the form. 


^■/=E n 

Wj} >^^=1 

j - 


MN „ 2 


MN 


^ 1 N H + ^*/ ) ’ (41) 


i=i 


where nj = n* — 1, n*, n* -f 1. We will see below that as the limit of iV —)■ oo is taken, the discrete 
modes in N will become continuous spectra in k over the Brillouin zone of the reciprocal space. 

Now we are ready to put all the phonon matrix elements together and perform the configura¬ 
tional sum. To do this we follow the steps of Huang and Rhys llllll . but generalize it for a system 
with multiple phonon frequencies. Eor multiple phonon bands, we assume that the frequency vari¬ 
ation within each band is much smaller than the frequency difference between the bands. This is 
the flat band approximation that is complemented with the requirement of finite spacing between 
the bands. We finally find, 

Fj = exp 


Ecothl 

Y 

J 

\ 1 

[ 2kT 

y2kT J 


smh.{tvjjj/2kT) 


(42) 


and 


F = 


Hi 


where 


tel I Vi=i / i=i 


J 

r 


g -'pj-i-i 

smhihojj/2kT' 


sinh(/ia;j/2A;T) 

S, 



^m\i{hi>jj/2kT) 



D{uj) 




Ej = l 

(43) 

(44) 


and Ip is the modified Bessel function of order p. 


B. Linear phonon matrix elements 


To evaluate the phonon matrix elements for the linear term, we rewrite it in terms of the normal 
mode coordinates Qj, 


E 

/ 


Y,M,{X,\q,\Xi) 


Y.Y.^M,{X,\qi\X,)\- 


f j 

1 V — 

2 ^ ^ 
/ 


(^/in [1 -f XMjQj exp{i(j)j)] \Xi) 


, (45) 


A=0 
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where (pj is a random phase introdueed to eancel out the cross terms, and, 




(46) 


The rest of the steps are exactly the same as for the zeroth order matrix elements. Using the 
integrals provided in Appendix UVl the matrix elements for the phonon part are. 


{Xn+i{qj + 6qj)\ [1 + AMjgjexp(i(/)j)] \Xn{qj)) = 


{n + l)a; 

w 


^ ^qj - -exp{i(j)j 


U-i 


{Xn-i{qj + Sqj) \ [1 + XMjqj exp{i(j)j)] |X„(gj)) = 


nuj 


6qj + 


XhMi 


Ui 


exp^icpj 


(47) 

(48) 


and. 


{Xn{qj + Sqj) \ [1 + XMjqj exp{i(j)j)] \Xn{qj)) = 1 - exp{i(j)j 


= 1 - 


2N 2 


XMj6qj exp(i0j). 


Define, 


(49) 


^±(A) = 


n + 1 


n 


2h^ 


N 


5qj =F 


XhMi 


bJn 


exp(i0j) 


n + 1 


UJ 


n 




exp 


XhM. , , , 

=F2—— exp(*(/)j) 


(50) 


The approximation in the second step is accurate to A^, with the ocnsideration that terms such as 
A^ sin 2(f)j and A^ cos 2(1)j drop out after the configurational average. Then, 


V';^ ^ + 1)5',-, 

^+(A) n + 1 
^_(A) n 

The A-dependent line-shape factor for a single phonon band is, 

c f H- 


XhMj , , . 
exY){i(i)j) 

UjSqj 


(51) 

(52) 


Fj{X) = exp 
exp 


-2Xp 


r — Sj coth 
hM. 


2kT """" V 2kT 


XN6qj I Mj exp {i(j)j ) | 






smh.{huj/2kT) 


X 


Ujdqj 


exp{i(f>j) 


(53) 


Let us now compare the two A factors by evaluating the ratio 

Nuj5q^ 


2h 


h 


(54) 
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For a hydrogenated vaeaney defeet our ealeulation shows that 6R 
Using m ^ 4.66 x 10“^® kg for the Si atom and uj ~ 10^^ seo“^, 


Nujjdq 


2h 


0.09. 


0.2 A for the nearest Si atom, 
we have, 


(55) 


Thus the first A faetor has a mueh smaller eontribution than the seeond one. The final linear phonon 
squared matrix element is, 


F = ^ V i — 

1 2Uk f- 1 dX^ 

{pjl 


M 

j=i 


A=0 




\ 1 

sinh(/ia;j/2/cT) 

smh.{hjjjj/2kT) 

ipj 


smh.{tvjjj/2kT) 


D(uj 


.(56) 




C. Ratio of zeroth-order and linear terms 


From the different expressions for the zeroth-order and the linear phonon matrix elements, we 
ean estimate the ratio between the linear term and the zeroth-order term in the transition rate. This 
is of the order of 


2 Mjhpj 

To estimate , we note that the leading term in Mj is (see Eq. (@1)), 


(57) 


0$ f 

Mj « ( 58 ) 

To estimate d^f/dqj, we assume rigid atomie orbitals, where the atomie wave funetions move 
rigidly in spaee with eaeh atom. The derivative of sueh a wave function with respect to atomic 
displacements simply reflects the change in the relative spatial phase, which is dictated by the 
phonon wave vector, 

$/exp(i0), (59) 

where \j is the acoustic wavelength for mode j, m is the mass of an atom, and 0 is the phase 
factor due to the movement of the atoms, which is different in each Born-von-Karman supercell. 
Integrating over all N Bom-von-Karman supercells, the sum of the exp(i0) factors scales as 1/N 
for large N. Thus, 
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Finally, pj is mostly zero, occasionally taking the values ±1, and 5qj ~ {m/N)6R where 6R is 

the largest atomic displacement and m is the mass of the corresponding atom. The ratio between 
the linear and zeroth order terms simplifies to, 

2 

(61) 


2 (—)■ 

V cm6R J 


where c is the sound velocity in the material. For a hydrogenated vacancy defect our calculation 
shows that 6R ~ 0.2 A for the nearest Si atom. Using this number and c ~ 8 x 10^ m/s for bulk 
silicon and m ~ 4.66 x 10“^® kg for the Si atom, we find. 


2 f ^ 3.6 X 10"^ 

\cmoR J 


(62) 


Thus the linear phonon term (non-adiabatic term) is several orders of magnitude smaller than the 
leading BOA term. 

D. Monte Carlo method for configurational sum 


The summation over all configurations {pj} involves a large number of terms when P = 
\Pj I greater than a few. We use a Monte Carlo approach to calculate this sum. For a given 
number of phonon modes, P, and a given number of bands, B, we use Monte-Carlo to construct 
a fixed number of configurations, K. We rewrite the sum over the configurations as a sum over 
the number of phonons P of a configuration, a sum over the number of bands B used to con¬ 
struct a configuration with P phonons and a sum over the configurations sampled (Monte Carlo 
steps). In each Monte Carlo step, we randomly pick B bands and then we construct all the possible 
configurations with P phonons constructed by these bands. 

In order to generate and count the configurations correctly, we first rewrite Eq. (|4^ as, 

, P K ( / M 


F = 




EE”®!]! 


X 


P=1 B=1 


{Pj}' 


\j=l 



^ Aj+i 

\ ] 

smh(hjjjj/2kT) 

smh.{hWj/2kT) 

S, 

smh{hUj/2kT) 




(63) 




Then, we normalize the sum so that the total weight, wp, in each sub-group of configurations 
(configurations with the same number of bands) is equal to the total number of possible configu- 
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rations for this number of bands, 


1 Ml 
K B\{M-B)\' 


(64) 


All configurations with up to four phonon modes are constructed and calculated explicitly. For 
configurations with more than four phonons, all the configurations constructed with up to three 
bands are calculated explicitly and the above equations are use to calculate the line shape function 
for configurations with more than three bands. 

The last step in the Monte Carlo scheme is to collect the line shape function into different en¬ 
ergy bins for a distribution. To do this, we note that with an incomplete sampling of the phase 
space via Monte Carlo, we may not be able to resolve the energy distribution to arbitrary accu¬ 
racy. Specifically, when we sample one configuration and weigh it according to Eq. (l64l) . we are 
effectively using it to approximate several configurations with different energies. Thus, the energy 
resolution must be consistent with the number of configuration samples - fewer configurations 
should correspond to coarser energy resolution. For this reason, we define the energy bin width 
separately for each value of P based on the requirement that there is at least one configuration 
inside each energy bin. To ensure the correct normalization, we rewrite the phonon density of 
states for band j as. 


D{uj) 


1 

'AE 


DWE = ll 


(65) 


where AE is the energy bin width and we assume that the phonon band is sufficiently flat so that 
it falls entirely within one energy bin. Then Eq. (l6^ becomes. 


E = 


1 

AE 


K 


M 


P=1B=1 {p,}/ L \i=l 


X 



^ ^Pj+^ 

\ 1 

si\ih.{hjjjj/2kT) 

smh.{hwj/2kT) 

J-Pj 

S, 

/2kT) 


( 66 ) 




The evaluation of the linear phonon terms is similar. 


V. APPLICATION TO A DEFECT IN SILICON 

In this paper, we will present only one application of the theory and computer codes for the 
capture cross section of a prototype defect in Si, namely a triply hydrogenated vacancy with a 
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bare dangling bond. Our purpose here is to demonstrate the feasibility of calculations, especially 
the first-ever calculation of the line-shape function that is converged with respect to the number 
of phonon modes that are used to construct random configurations whose energy is equal to the 
amount of energy that needs to be dissipated following the instantaneous electronic transition. We 
defer calculations for defects for which experimental data are available to a future paper where we 
anticipate using hybrid functionals in the DFT calculations of the electronic matrix elements. Such 
calculations are computationally demanding, but would provide more accurate transition energies 
and electronic matrix elements. In addition, we plan to code the additional contributions from the 
linear terms which we estimated to be significantly smaller because they scale with the inverse of 
the mass of a typical atom in the defect cluster. It will be interesting to see how the two terms in 
the square brackets in eq. [T^add or subtract for different defects. 

In Fig. I, we show the values of calculated electronic matrix elements as a function of energy. 
At each energy value, there are a number of k points that contribute. Their contributions are 
indicated by red symbols. The size of the energy bin is determined by the number of k points. 
For the example shown in Fig. 1 the average matrix element as a function of energy is shown 
by the blue line. The size of the energy bin fixes the resolution. A smooth curve can only be 
obtained with very small energy bins, which requires a very large number of k points. It is clear 
from the figure that the capture electronic matrix element is relatively constant as a function of 
energy, whereby it seems best at this point to take it to be a constant, either an average value or the 
value at the threshold for capture, which introduces an error bar of a factor of ~ 1.7 (clearly, to 
validate the theory against accurate experimental data, we need a very accurate calculation in the 
near-threshold region). 

In Fig. 2, we show the calculated capture cross section using a constant matrix element to show 
clearly the convergence of the line-shape function as we increase the number of phonon modes 
that are used to construct configurations (the electronic matrix element is just a multiplier that sets 
the absolute value). The dominant contribution to the line-shape function comes from the balance 
between the modes with largest general coordinate displacement (GCD) and the growth of the 
number of allowed combinations with smaller GCD. Note that the curves are smooth because we 
employ millions of configurations at each energy and therefore we have very tiny energy bins. It is 
clear that a single-phonon-mode approximation would be very poor indeed. In Fig. 3 we show the 
convergence of the capture cross section at threshold (for electrons at the bottom of the conduction 
band), which is what is usually measured. Once more, it is clear that the single-phonon mode 
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Electron energy (eV) 


FIG. 1: Calculated electronic matrix elements as a function of the initial state electron energy for a triply 
hydrogenated vacancy in Si with a bare dangling bond. Red points: matrix element values at each energy 
for different k points; blue curve: Averaged matrix element over all k points for each energy. 


approximation would be inadequate. 

For a calculation of the cross section using electronic matrix elements that depend on energy, 
the resolution is limited by the energy bin size. We show the result in Fig. 4. Clearly, the size of 
the energy bin is important. For capture cross sections, one is often interested only in the thresh¬ 
old value. The calculations presented here are a prelude to calculations of hot-electron inelastic 
multiphonon scattering, for which the energy dependence is important. The energy dependence is 
also important in luminescence curves, Le., the classic Huang-Rhys problem that was treated in 
the single-phonon approximation in Ref. l23| (in the case of luminescence, MPPs dissipate only the 
relaxation energy of the defect, when one expects the phonon mode corresponding to the actual re¬ 
laxation to dominate; nevertheless, a fully convergent calculation would be needed to establish the 
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FIG. 2: Calculated electron capture cross section using a constant electron matrix element and different 
number of phonon modes. 

degree of aeeuraey one obtains with the single-mode approximation). The aeeuraey of the ealeu- 
lation of the line shape funetion is eontrolled by the aeeuraey of the ealeulation of the generalized 
displaeements. The latter depends on the aeeuraey of the calculation of the atomic displacements. 
We found that accuracy is enhanced significantly if we allow the entire supercell to relax, which 
allows the defect’s neighbors to relax more freely. At the same time, a dense fc-point mesh is nec¬ 
essary. In Fig. 5, we present the atomic displacements of the triply-hydrogenated Si vacancy as a 
function of the distance from the vacancy site for a 64-atom supercell. Using only one fc-point and 
not allowing the supercell to relax we get only the Si-atom near the defect to move significantly 
while the rest of the crystal remains essentially frozen (blue dots). This kind of relaxation leads to 
only a few phonon modes being significant and thus the system is artificially able to dissipate en¬ 
ergy efficiently at certain frequencies. On the other hand the well-relaxed crystal of the (3x3x3) 


23 












4 5 6 7 8 9 10 11 12 

Number of phonon modes 


FIG. 3: Convergence of the calculated electron capture cross section at the threshold as a function of the 
number of phonon modes. 

/c-points grid (red dots) has more atoms contributing to the generalized displacements and thus 
almost all the phonon modes contribute in the dissipation to the energy of the incoming electron. 
The use of supercells with more than 64 atoms would be prohitively expensive for the line-shape 
function calculation. 

VI. SUMMARY 

We have presented a comprehesive theory of inelastic multiphonon carrier capture and scatter¬ 
ing processes. We showed that, under non-equilibrium conditions, i.e., in the presence of currents 
or hot electrons, the defect potential is primarily responsible for capture throught a zeroth-order 
term in an expansion in terms of the atomic displacements (relaxation) that accompanies capture. 
These terms were not included in any prior theory. Instead, the focus has always been on the linear 
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Electron energy (eV) 


FIG. 4: Calculated full capture cross section using the electron matrix element from Fig. 1 and 12 phonon 
modes in the line-shape function. 

terms, which we showed here to be much smaller because they depend on the inverse of the mass 
of typical atoms in the defect complex. The linear terms are dominant only in the limit of ther¬ 
mal equilibrium. For the first time, we used accurate all-electron wave functions obtained by the 
PAW method for the electronic matrix elements and an accurate Monte Carlo scheme to sample 
random configurations of up to 12 distinct phonon modes for the line-shape functions to achieve 
convergence (a single-phonon-mode approximation has been standard in prior calculations). We 
presented results for a prototype defect. More accurate hybrid exchange-correlation functionals 
are needed to produce results that are accurate enough for comparison with experimental data. In 
addition, a reliable comparison with data can only be made with experimental measurements of 
capture cross sections simultaneously with the determination of the elastic mean-free-path and the 
capture mean-free-path, as they appear in Eq. ([T])- 
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FIG. 5: Atomic displacements of the triply-hydrogenated Si vacancy as a function of the distance from the 
vacancy site for a 64-atom supercell. 
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I. EVALUATION OF THE ELECTRONIC MATRIX ELEMENT FOR BOA TRANSITION 


In the basis of |^n)> the unperturbed Hamiltonian is diagonal with eigenenergies e„. The 
total electron Hamiltonian H = has coupling terms only between states |\hi) and |'k/). 

We can, therefore, construct solutions of 


+ = E|$). 


( 1 ) 


in the form |$) = al'l/j) + so that 


AH 


AH, 

fi 


if 


a 

b 


= E 


( 2 ) 


There are two sets of solutions. 


^i{f) - 


Ci + 6/ ± (cj - e/)2 + A\AHij\ 


( 3 ) 


where state i takes the + sign and state / take the — sign, since Ei > Er. The coefficients satisfy 


€i(ii T AHijbi Ej^Oji 


and 


— 1- 

There is an arbitrary phase factor within Oj. We can define a set of solution as. 


( 4 ) 


( 5 ) 


ai = b*f = 




AH 


if 


Ei — Ef 


and 




Ei-E 


f 


( 6 ) 


( 7 ) 


'iHp 


AH, 


if 


\Ei-Ef\ 

If we can compute the overlap integral ($/|'kj) = af, then we can solve for | Aifj/p from |a/p 


and find. 






2'i2 if' 


( 8 ) 
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To be consistent with the phase of Eq. (|7]), we have, 



( 9 ) 


The wave function |\I/i) is related to that of a perfect crystal I'P®) through a perturbation expan 
sion, 



( 10 ) 


Because Hi has only nonzero elements between the states |^i) and |^/), for j ^ i, f, the wave 
functions |^j) = |$j) so that = 0. Thus, to first order in the defect potential. 




( 11 ) 


and, assuming that |($/|\1'°)| -C 1, we arrive at Eq. (|2^ . which simplifies the evaluation of the 
overlap integral. 

II. EVALUATION OF THE GRADIENT TERMS 

Using the result in the previous section for the matrix element we now calculate the 

gradient terms in Eq. (fT9l) . f\Vi) ■ Neglecting higher order 

terms, the first gradient term is. 


= - Gif - f), ( 1 ) 


where in the last step we used the fact that VRj.ej = 0 (the initial state is at equilibrium) and the 
Helmann-Eeynman theorem for VRj.e/. Erom Eq. (fT^ we have. 



( 2 ) 


where we used VR^^ifei = VR^^ife- Because lihi/) = for i' ^ i, f and = 1 + 

0(|($/|^'i)P), we have. 







(3) 


Similarly, 





(4) 
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Combining these results and noting that does not have diagonal components, we arrive at 

VR.MfO + = {4,|VR.i/|4i> - (5) 

We can use Eq. (fTTI) and approximate |^/) ^ |<h/) to get Eq. (l24l) . 


III. EVALUATION OF THE OVERLAP INTEGRAL WITHIN THE PAW 

Consider the problem of evaluating the overlap integral between two wave functions 

from two different solids (e.g., one is a perfect crystal and the other contains a defect). Using the 
PAW expansion of the full wave functions: 

= + ( 1 ) 

where |\^) is the pseudo wave function and |\l'^'®)a and are the atomic wave functions 

inside the augmentation sphere of each atom a, and similarly, 

1$) = ||.) + (2) 

Now, (^'1$) is given as: 

(^ 1 $) = (^(^1 j ^ 

The first term, (^1$), is the overlap of the pseudo wavefunctions and can be easily calculated 
since the pseudo wavefuntions are expanded in the same base set of plane waves. 

In order to evaluate the terms and —a we make use 

of the unitary operators constructed by the projectors \p) and the pseudo atomic wavefunctions 

l«: 

b,ib 

and 

= l (5) 


29 


inside the augmentation sphere of eaeh atom b of the perfect crystal and each atom a of the solid 
with the defect respectively. Thus: 

b,ib 

and 

a,ia 

Equations ® and (|7]) ensure that in the case that if the two solids are identical, i.e. |^) and 
1$) are eigenstates of the same Hamiltonian and the augmentations spheres are identical, the one 
center expansion |0)(p|^) of tho pseudo wavefunction is identical to the pseudo wavefunction 
1^) inside the augmentations sphere and 

( 8 ) 

To evaluate Eqs. ® and ([7]), we need the projections of the pseudo wavefunctions of the first 
solid to the projectors of the atomic wavefunctions of the second solid, and vise versa for 

the projections {pf^ | $). This can be easily calculated since both the pseudo wavefunctions and the 
projectors are expanded in the same base set of plane waves. 

The difficulty in evaluating the last term in Eq. ([3]) —a (^^‘^1) is 

that the cutoff spheres for the two wave functions are usually not identical. We can bypass this 
difficulty by evaluating the integral with the assistance of a complete set of plane waves |k), 

-a (^^^1) - \^f’^)b) = 

k 

= (^'^''Ik)) ((k|$^®)6 - (k|<h^^),) . 

k 

The plane waves can be expanded in either sphere as 

= 4^5^ i‘Mkr)Y]:,Xk)Y,Ur). (9) 

Im 

and using, 

|k> = (10) 

the all-electron and the pseudo atomic wave functions is written as: 

|4“>, = 5^ iji.). (11) 
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.(>P-"®|k) {>l.™|k) = ^(>l'|fi„.>e'‘‘'«-i'“y,:„.(k) / 3 iAkT)(Ril - 

V t' a.i. ■'» 

(13) 


Att ~ - r’’ 

(k|4-“)t - {k|4f®)t = ^ 3,,(kr)(Rff^ - R^fydr, 


IV. PHONON INTEGRALS 


The overlap matrix between the initial and final states for the mode j is, 

T ^Qj}\^nk{Qj)) j ^R-j+Pji^j ^Qj}^nj{Qj}dqj. (1) 

where n* = Uj and n| = Uj + Pj. 

For convenience, we drop the subscript j for rij and pj. Expanding Xn{qj + 5qj) in terms of 

5qj, 






Defining the raising and lowering operators 


we have. 


h d [uJj 

“* “ ’'V ^ V 2fi*’ 


5'+^n(*?j) \/Ti “h ^ 




Subtracting the two, we find, 


^Xn{qj) = \[^{d- - d+)Xn{qj) = ^[^Xn-i{qj] 


\n + l)ujj 


Xn+i{qj). (6) 


Using this recursive relation, we find that the lowest order term for J Xn{qj + 6qj)Xn+k{qj)dqj is 
dqf^. Therefore, for small 6qj only k = ±1 terms dominates. It means that each mode would at 
most emit or absorb a single phonon. 
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The result for the integrals are, 




dq 


2h 


(note that this was ineorreetly given as —{^yhmu/2)y/n + 1 in Ref. [nl) . and, 

f dXn{q) 


dq 


-Xn-i{q)dq = 


nuj 


For linear phonon matrix elements, we have. 


qjXniqj) — W T —(d- + d+)Xn{qj) — J- — Xn-i{qj) + \ — ^—^n+iidj] 


2uj4 


The integrals needed are. 


2uji 


2oj4 


Xn{q)Xn+i{q)qdq = 


{n + l)h 


2uj4 


( 7 ) 


( 8 ) 


( 9 ) 


( 10 ) 


and. 


Xn{q)Xn-i{q)qdq = 


nh 

2u}j 


dXn{q) V ( \ ^ ^ 

—^Xn{q)qdq = 


( 11 ) 


( 12 ) 


V, LINE SHAPE FUNCTION 

We first eonsider a single phonon band, i.e., all of the phonon modes ujj = Ci;(kj) form a 
single eontinuous band deseribed by wave veetors kj. Beeause of the Born-von-Karman periodie 
boundary eondition, the phonon band is diseretized into N modes. Suppose that s modes go down 
by one quantum and s + p modes go up by one quantum. Then the line shape faetor, Eq. (|4TI) with 
M = 1, eontains eontributions formed from the following produets, 
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where X]j=i ~ ^j) i® the energy differenee beeause of the different phonon frequeneies of 

the initial and final eonfiguration of the defect and tj, /_, and /+ are defined as: 
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( 2 ) 


Xni{qi)Xni{qi + ^Ql)d<ll 


A naive way to sum over all possible configurations is to neglect the difference in the frequencies 
and apply the same counting method as Huang and Rhys 111 ill] to write the configurational sum for 
all such combinations of phonons as, 
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(3) 


This would not be correct if the frequencies are different for each mode. Furthermore, the sum¬ 
mation over configurations for large N is needed to integrate out the 6 function. Therefore the 6 
function cannot be left outside the summations. Let us consider one term in the 5 function at a 
time. Consider one the plus terms hiof^ and insert the 6 function into one of the summations. 
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For large N, each of the summations inside the curly brackets can be converted into integrals and 
evaluated, 

n + 1 
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where f2k is the volume of the reeiprocal space Brillouin zone. In the last step we assumed that 
the frequency and displacement do not change with k. 

In order to evaluated the last factor that includes the S function, we note that each term in the 
summation over m has a different which spans the entire phonon band when m scans from 1 
to N. Thus as we convert the sum over m to integral over k, the argument is also converted to 

Wk, 


N 
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<{ [^L+ Y W - 5^ 
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^ [ /k,+5 (f^l+ Y 
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ojj - oOj) + eif dk 


= S. 


D{u) 




( 6 ) 


^■''+Ei6s+p-i^i ‘^i)+ei/=o 


where D{oj) is the phonon density of states. Combining the above equations and then setting all 
frequencies to ce, Eq. dH) now becomes, 


N 


s!(s +p)\ 1 

K 7 = 1 


JJqsiS‘+’’ 


DH 




(7) 


P^+EEi n'.hiujf. -uj’j)+eif=0 


But there is one such contribution for each Uk or oji in the 5 function, regardless of the sign of the 
frequency. For s modes subtracting a phonon and s + p modes adding a phonon there are total 
2s + p such contributions. We thus sum over all the terms and obtain. 
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Finally, the factor 11^1 i®’ 
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n 

i=i 


1 - 


{2n + l)a; ^ 2 


H 2N 


Ah 


-6q^ 


= exp [-(S'++ S'-)]. (9) 


The line shape factor for a single phonon band is. 
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To generalize the above expression to multiple phonon bands, the normalization faetor must be 
evaluated with a summation over both the band index and the k points within eaeh band. If we use 
Fj to denote the faetor for a band that adds net pj phonons, i.e.. 
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( 11 ) 


then in a similar manner as for the ease of a single phonon band, Fj is evaluated to be. 
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—j exp[-Sj{2nj + l)]Ipj 
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( 12 ) 


Now we insert the 6 funetion into the produet of Fj in the same manner as in the ease of a 
single band to form the full line shape faetor, one phonon band at a time. For now let us eonsider 
the ease where all p/s are positive. We have. 
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(13) 


where j" is one of the phonon bands and we have used Eqs. dS]) and db]). Summing over all possible 
j" terms and with an additional summation over all eonfigurations {pj}, we find. 
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If some of the p/s are negative, we need to switeh the roles of 5*+ and S', following Ref. ill. 
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Redefining Sj + pj Sj and Sj Sj — pj in Eq. (fTSl) . the factor corresponding to pj becomes, 
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(15) 


using the recurrence relation for the Bessel functions. Therefore Eq. (fT4l) is valid for both positive 
and negative p/s. Applying thermodynamic average to the occupation numbers, rij is replaced by 
the Bose-Einstein distribution function. 
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and 


we obtain Eqs. (1421) and 
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